##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
suppressWarnings(library(edgeR, quietly = T))
library(RColorBrewer)
library(ggplot2)
library(ggrepel)
##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--out_path"), type = "character"),
    make_option(c("--foldchange_t"), type = "character"),
    make_option(c("--q_t"), type = "character")
)

if(1!=1){
    
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/RSEM/CombineCounts.FilterLowExpression-MergeMutiSample.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene"
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file
q_t <- as.numeric(opt$q_t)
foldchange_t <- as.numeric(opt$foldchange_t)

dir.create(out_path , recursive = T)

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_tpm <- data.frame(fread(rsem_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")

class_type <- c( "GS" , "CIN" , "MSI")
##########################################################################################

count_norm <- dat_tpm
rownames(count_norm) <- count_norm$gene_id
count_norm <- count_norm[,-which(colnames(count_norm) == "gene_id")]
##只提取IM的样本
count_norm <- count_norm[,grep( "_DGC", colnames(count_norm) , value = T)]
colnames(count_norm) <- paste0(sapply(strsplit(colnames(count_norm),"_"),"[",1),"_",info$Molecular.subtype[match(sapply(strsplit(colnames(count_norm), "_"), "[", 1), info$Tumor)])
print(colnames(count_norm))

conditions <- factor(sapply(strsplit(colnames(count_norm) , "_") , "[" , 2))
table(conditions)
##########################################################################################
## 多组，两两做差异表达
result_diff <- c()

for( k in 1:(length(class_type)-1) ){
    class1 <- class_type[k]
    print(class1)

    for( j in (k+1):length(class_type) ){
        class2 <- class_type[j]
        print(class2)

        ## Run the Wilcoxon rank-sum test for each gene
        use_sample <- grep( paste0(class1 , "|" , class2) , colnames(count_norm) , value = T)
        count_norm_use <- count_norm[,use_sample]
        conditions_use <- grep( paste0(class1 , "|" , class2) , conditions , value = T)

        pvalues <- sapply(1:nrow(count_norm_use),function(i){
            data<-cbind.data.frame(gene=as.numeric(t(count_norm_use[i,])),conditions_use)
            p <- wilcox.test(gene~conditions_use, data)$p.value
            return(p)
        })
        fdr=p.adjust(pvalues,method = "fdr")

        ## Calculate the fold-change for each gene
        conditionsLevel <- c(class1 , class2)
        dataCon1 <- data.frame(count_norm_use[,c(which(conditions_use==conditionsLevel[1]))])
        dataCon2 <- data.frame(count_norm_use[,c(which(conditions_use==conditionsLevel[2]))])
        print(head(dataCon2))
        if(ncol(dataCon1) > 1) {
            row_means1 <- rowMeans(dataCon1)
            } else if(ncol(dataCon1) == 1) {
                row_means1 <- dataCon1[[1]]
            }
        if(ncol(dataCon2) > 1) {
            row_means2 <- rowMeans(dataCon2)
            } else if(ncol(dataCon2) == 1) {
                row_means2 <- dataCon2[[1]]
            }
        foldChanges <- log2(row_means2/row_means1)

        ## Output results based on FDR threshold
        outRst<-data.frame(log2foldChange=foldChanges, pValues=pvalues, FDR=fdr)
        rownames(outRst)=rownames(count_norm_use)
        outRst=na.omit(outRst)

        ## class1代表分子
        outRst$class1 <- class2
        outRst$class2 <- class1
        outRst$meanTMM_class1 <- rowMeans(dataCon2)
        outRst$meanTMM_class2 <- rowMeans(dataCon1)
        outRst$gene_id <- rownames(outRst)

        result_diff <- rbind( result_diff , outRst)
    }
}

##########################################################################################

colnames(result_diff) <- c("log2FoldChange" , "pvalue" , "padj" , "class1" , "class2" , "meanTMM_class1" , "meanTMM_class2" , "gene_id")

image_name <- paste0( out_path , "/DiffGene_DGC.tsv" )
write.table( result_diff , image_name , row.names = F , sep = "\t" , quote = F )

##########################################################################################
dat_diff <- result_diff
dat_diff <- merge( dat_diff , dat_gtf , by = "gene_id" )
##########################################################################################
diffplot <- function(use_dat = use_dat , image_name = image_name,lauren,Class1,Class2){

    #对原数据进行处理
    use_dat$log10_q_Value <- -log10(use_dat$padj)
    use_dat$log_FC <- use_dat$log2FoldChange
    use_dat$gene <- use_dat$Hugo_Symbol

    data <- use_dat[,c('gene','log_FC','log10_q_Value')]
    colnames(data) <- c('gene','log_FC','-log10_P_Value')

    #设置阈值
    logFC_cutoff <- log2(foldchange_t)
    log10_P_Value_cutoff <- -log10(q_t)

    options(ggrepel.max.overlaps = 20)

    plot1 <- ggplot(data = data,aes(x = log_FC,y = `-log10_P_Value`))+
      geom_point(data = subset(data,abs(log_FC)<logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'gray',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)<log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'gray',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC>logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'red',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC< -logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'darkgreen',alpha = 0.4)+
      theme_bw()+
      theme(legend.title = element_blank(),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position = 'none',
            axis.line = element_line(colour = "black"))+
      labs(x='log2(fold change)',y='-log10(adjusted p-value)')+
      geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),lty = 3,col = 'black',lwd = 0.4)+
      geom_hline(yintercept = log10_P_Value_cutoff,lty = 3,col = 'black',lwd = 0.4) +
      geom_text_repel(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                            aes(label = gene),size = 3,col = 'black')+
      ggtitle(paste0(lauren," : ",Class1," vs ",Class2))+
       theme(plot.title = element_text(size = 16))
    ggsave( image_name , plot1 , width = 10 )
}


use_dat <- subset( dat_diff , class2=="GS" & class1=="CIN")
image_name <- paste0( out_path , "/DGC_CIN_GS.DiffGene.pdf" )
diffplot(use_dat = use_dat , image_name = image_name,lauren="DGC",Class2="GS" , Class1="CIN")

use_dat <- subset( dat_diff , class2=="CIN" & class1=="MSI")
image_name <- paste0( out_path , "/DGC_MSI_CIN.DiffGene.pdf" )
diffplot(use_dat = use_dat , image_name = image_name,lauren="DGC",Class2="CIN" , Class1="MSI")

use_dat <- subset( dat_diff , class2=="GS" & class1=="MSI")
image_name <- paste0( out_path , "/DGC_MSI_GS.DiffGene.pdf" )
diffplot(use_dat = use_dat , image_name = image_name,lauren="DGC",Class2="CIN" , Class1="MSI")